# Generate Figure 3 in the paper 
# set seed
set.seed(102112516) 

N <- 15000 # sample size of pseudo sample
m <- N     # sample size of pseudo sample

# lognormal parameters
s1 <- 1
s2 <- 2
p1 <- 0.2
p2 <- 0.8


X <- vector(mode = "list", length = 6)

# Generate samples
X[[1]] <- exp(mvrnorm(N,c(0,0),matrix(c(s1,p1,p1,s1),nrow=2)))
X[[2]] <- exp(mvrnorm(N,c(0,0),matrix(c(s1,p2,p2,s1),nrow=2)))
X[[3]] <- exp(mvrnorm(N,c(0,0),matrix(c(s2,p1,p1,s2),nrow=2)))
X[[4]] <- exp(mvrnorm(N,c(0,0),matrix(c(s2,p2,p2,s2),nrow=2)))
X[[5]] <- exp(mvrnorm(N,c(0,0),matrix(c(s1,p1,p1,s2),nrow=2)))
X[[6]] <- exp(mvrnorm(N,c(0,0),matrix(c(s1,p2,p2,s2),nrow=2)))

# Obtain vector quantiles
ot1 <- vquantile(X[[1]],rep(1/N,N))
ot2 <- vquantile(X[[2]],rep(1/N,N))
ot3 <- vquantile(X[[3]],rep(1/N,N))
ot4 <- vquantile(X[[4]],rep(1/N,N))
ot5 <- vquantile(X[[5]],rep(1/N,N))
ot6 <- vquantile(X[[6]],rep(1/N,N))

# Obtain ILF to find alpha-Lorenz curves
ilf1 <- ILF(X[[1]],rep(1/N,N),ot1,m)
ilf2 <- ILF(X[[2]],rep(1/N,N),ot2,m)
ilf3 <- ILF(X[[3]],rep(1/N,N),ot3,m)
ilf4 <- ILF(X[[4]],rep(1/N,N),ot4,m)
ilf5 <- ILF(X[[5]],rep(1/N,N),ot5,m)
ilf6 <- ILF(X[[6]],rep(1/N,N),ot6,m)

# Generate a grid
sequence <- rep(seq(from = 0, to = 1, by = 1/(nrow(ilf1$estimate) - 1)))
s <- 0
t <- 0
x <- rep(0,length(sequence)^2)
y <- x
for (i in 1:length(sequence)){
  s <- s+1
  k <- length(sequence)*(s - 1) + 1
  j <- length(sequence)*s
  x[k:j] <- rep(sequence[s],length(sequence))
  y[k:j] <- sequence
}

# Calculate the coordinates of the corners of the identical distribution
# There are more alphas than necessary, we only use 0.9 in the figure.
alpha <- c(0.1,0.25,0.3,0.4,0.5,0.6,0.7,0.8,0.9)
egal <- rep(0,9)
for (k in 1:9){
  # This formula is given in paper section 1.2.2 example 1, set d = 2
  egal[k] <- uniroot(function(x){x*(1-log(x)) - alpha[k]},c(0.001,1))$root
}

lr <- data.frame(x,y,c(ilf1$estimate),
                 c(ilf2$estimate),
                 c(ilf3$estimate),
                 c(ilf4$estimate),
                 c(t(ilf5$estimate)),
                 c(t(ilf6$estimate))
)

names(lr) <- c("x","y","l1","l2","l3","l4","l5","l6")

levels <- 0.9
sz <- 1.25

lcontours <- ggplot(lr,aes(x=x,y=y)) + 
  annotate("text", x=0.65, y=0.15, label= expression(italic(z)["1"]), family = 'serif',size = 7) +
  annotate("text", x=0.15, y=0.65, label= expression(italic(z)["2"]), family = 'serif',size = 7,angle = 270) +
  geom_contour(aes(z=l1,color = "1989",linetype = "1989"),breaks = levels, size =sz) + 
  geom_contour(aes(z=l2,color = "1992",linetype = "1992"),breaks = levels, size = sz) + 
  geom_contour(aes(z=l3,color = "1995",linetype = "1995"),breaks = levels, size = sz) + 
  geom_contour(aes(z=l4,color = "1998",linetype = "1998"),breaks = levels, size = sz) + 
  geom_contour(aes(z=l5,color = "2001",linetype = "2001"),breaks = levels, size = sz) + 
  geom_contour(aes(z=l6,color = "2004",linetype = "2004"),breaks = levels, size = sz) + 
  annotate("segment",x = 0, y = 0, xend = 0.6, yend = 0.6,linetype = 3, size=1)+
  scale_color_manual(
    name = expression(paste("(",sigma[1]^2,", ",sigma[2]^2,", ",rho,")")),
    labels = c("(1,1,0.2)",
             "(1,1,0.8)",
             "(2,2,0.2)",
             "(2,2,0.8)",
             "(1,2,0.2)",
             "(1,2,0.8)"),
    values=c("#ee99aa",
             "#994455",
             "#6699cc",
             "#004488",
             "#eecc66",
             "#997700")) +
  scale_linetype_manual(
    name = expression(paste("(",sigma[1]^2,", ",sigma[2]^2,", ",rho,")")),
    labels = c("(1,1,0.2)",
               "(1,1,0.8)",
               "(2,2,0.2)",
               "(2,2,0.8)",
               "(1,2,0.2)",
               "(1,2,0.8)"),
    values=c("solid",
             "solid",
             "dashed",
             "dashed",
             "dotdash",
            "dotdash")) +
  coord_cartesian(xlim = c(0.15, 0.65),ylim = c(0.15, 0.65)) +
  theme(plot.title = element_text(hjust = 0.5)) + xlab("")+ ylab("")+
  theme(panel.border = element_rect(colour = "black", fill=NA, size=0.6),
        panel.background = element_rect(fill = 'transparent'),
        panel.grid = element_blank(),
        legend.position=c(0.76,0.76),
        legend.key = element_rect(fill = 'transparent'),
        text = element_text(family = "serif",size=22))+
  scale_y_continuous(breaks = seq(0.2,1,0.1),position = "right",labels = seq(0.2,1,0.1)) +
  scale_x_continuous(breaks = seq(0.2,1,0.1),position = "top",labels = seq(0.2,1,0.1)) +
  theme(legend.key.width=unit(3,"line"))
ggsave(plot = lcontours, width = 6, height = 6, filename = "fig3.png")






